#include "funcRichards.h"

using namespace QPCRModel;

funcRichards::funcRichards(){
    a = 1;
    b = 1;
    x0 = 1;
    y0 = 0;
    d = 1;
}

int funcRichards::getParamNum(){
    return 5;
}

double funcRichards::f(double x){
    return y0 +  a / pow(1 + pow(M_E, (x0-x)/b), d);
}

double funcRichards::g(double y){
    if(y < y0 || y > y0 + a){
        return 0;
    }
    return x0 - b * log( pow(a/(y-y0), 1/d) -1);
}

double funcRichards::loss(double x, double y){
    return y - ( y0 + a / pow( 1 + exp((x0-x)/b), d ) );
}

double funcRichards::loss(double x, double y, const Eigen::VectorXd &X){
    return y - ( X[3] + X[0] / pow( 1 + exp((X[2]-x)/X[1]), X[4] ) );
}

double funcRichards::loss2(double x, double y){
    return pow(y - ( y0 + a / pow( 1 + exp((x0-x)/b), d ) ), 2);
}

double funcRichards::loss2(double x, double y, const Eigen::VectorXd &X){
    return pow(y - ( X[3] + X[0] / pow( 1 + exp((X[2]-x)/X[1]), X[4] ) ), 2);
}

double funcRichards::diff1(double x){
    return (a*d*exp(-(x - x0)/b))/(b*pow(exp(-(x - x0)/b) + 1, d+1));

    //return (a*d*exp(-(2*(x - x0))/b)*(d + 1))/(b*b*pow(exp(-(x - x0)/b) + 1, d+2)) - (a*d*exp(-(x - x0)/b))/(b*b*pow(exp(-(x - x0)/b) + 1, d+1));
}

double funcRichards::diff2(double x){
    return (a*d*exp(-(2*(x - x0))/b)*(d + 1))/(b*b*pow(exp(-(x - x0)/b) + 1, d+2)) - (a*d*exp(-(x - x0)/b))/(b*b*pow(exp(-(x - x0)/b) + 1, d+1));
}

Eigen::VectorXd funcRichards::jacobian(double x, double y){
    Eigen::VectorXd res(5);
    res <<  -1/pow(exp((x0 - x)/b) + 1, d),
            -(a*d*exp((x0 - x)/b)*(x0 - x))/(b*b*pow(exp((x0 - x)/b) + 1, d+1)),
            (a*d*exp((x0 - x)/b))/(b*pow(exp((x0 - x)/b) + 1, d+1)),
            -1,
            (a*log(exp((x0 - x)/b) + 1))/pow(exp((x0 - x)/b) + 1, d);
    return res;
}

Eigen::VectorXd funcRichards::jacobian2(double x, double y){
    Eigen::VectorXd res(5);
    // 好像用不上....不算了
    return res;
}


void funcRichards::setV(const Eigen::VectorXd &v){
    a = v[0];
    b = v[1];
    x0 = v[2];
    y0 = v[3];
    d = v[4];
}

Eigen::VectorXd funcRichards::getV(){
    Eigen::VectorXd res(5);
    res << a, b, x0, y0, d;
    return res;
}

double funcRichards::getLinearStageFunc(double *lk, double *lb){
    *lk = (a+y0) / b * pow(d/(d+1), d+1);
    *lb = -(*lk) * (x0 + b*log(d)) + (a+y0) * pow(d/(d+1), d) + y0;
    return *lk;
}


double funcRichards::curvity(double x, double norm_coef){
    if(norm_coef <= 0){
        norm_coef = a;
    }
    double res = -((norm_coef*d*exp(-(x - x0)/b))/(b*b*pow(exp(-(x - x0)/b) + 1, d+1)) - (norm_coef*d*exp(-(2*(x - x0))/b)*(d + 1))/(b*b*pow(exp(-(x - x0)/b) + 1, d+2)))/pow((pow(norm_coef,2)*d*d*exp(-(2*x - 2*x0)/b))/(b*b*pow(exp(-(x - x0)/b) + 1, 2*d + 2)) + 1, 1.5);
    //double res = -((a*d*exp(-(x - x0)/b))/b*b*pow(exp(-(x - x0)/b) + 1, d + 1) - (a*d*exp(-(2*(x - x0))/b)*(d + 1))/(b*b*pow(exp(-(x - x0)/b) + 1, d + 2)))/pow((a*a*d*d*exp(-(2*x - 2*x0)/b))/(b*b*pow(exp(-(x - x0)/b) + 1,2*d + 2)) + 1, 3/2);
    return res;
}

std::pair<double, double> funcRichards::getRange(int i){
    if(i == 0){
        return std::pair<double, double>(0, INFINITY);
    }else if(i == 1){
        return std::pair<double, double>(0.3, 10);
    }else if(i == 2){
        return std::pair<double, double>(0, 60);
    }else if(i == 3){
        return std::pair<double, double>(-1000, INFINITY);
    }else if(i == 4){
        return std::pair<double, double>(0.1, 1000);
    }else{
        return std::pair<double, double>(-INFINITY, INFINITY);
    }
}


void funcRichards::init(double SPE, double EPE, int cycleNum, double ymin, double ymax){


    // 根据findExp2函数，SPE和EPE定义为线性期两侧斜率小于2的点
    if(EPE > cycleNum){
        a = (ymax - ymin) * 2;
    }else{
        a = (ymax - ymin);
    }

    b = 2;
    x0 = (SPE + EPE) / 2;
    y0 = ymin;

    d = 1;

}
